home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / fft / c_radix2.c < prev    next >
Encoding:
Text File  |  2001-05-14  |  7.6 KB  |  317 lines
  1. /* fft/c_radix2.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. int
  21. FUNCTION(gsl_fft_complex,radix2_forward) (TYPE(gsl_complex_packed_array) data,
  22.                       const size_t stride, const size_t n)
  23. {
  24.   gsl_fft_direction sign = forward;
  25.   int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
  26.   return status;
  27. }
  28.  
  29. int
  30. FUNCTION(gsl_fft_complex,radix2_backward) (TYPE(gsl_complex_packed_array) data,
  31.                        const size_t stride, const size_t n)
  32. {
  33.   gsl_fft_direction sign = backward;
  34.   int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
  35.   return status;
  36. }
  37.  
  38. int
  39. FUNCTION(gsl_fft_complex,radix2_inverse) (TYPE(gsl_complex_packed_array) data,
  40.                       const size_t stride, const size_t n)
  41. {
  42.   gsl_fft_direction sign = backward;
  43.   int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
  44.  
  45.   if (status)
  46.     {
  47.       return status;
  48.     }
  49.  
  50.   /* normalize inverse fft with 1/n */
  51.  
  52.   {
  53.     const ATOMIC norm = 1.0 / n;
  54.     size_t i;
  55.     for (i = 0; i < n; i++)
  56.       {
  57.     REAL(data,stride,i) *= norm;
  58.     IMAG(data,stride,i) *= norm;
  59.       }
  60.   }
  61.  
  62.   return status;
  63. }
  64.  
  65.  
  66.  
  67. int
  68. FUNCTION(gsl_fft_complex,radix2_transform) (TYPE(gsl_complex_packed_array) data,
  69.                         const size_t stride, 
  70.                         const size_t n,
  71.                         const gsl_fft_direction sign)
  72. {
  73.   int result ;
  74.   size_t dual;
  75.   size_t bit; 
  76.   size_t logn = 0;
  77.   int status;
  78.  
  79.   if (n == 1) /* identity operation */
  80.     {
  81.       return 0 ;
  82.     }
  83.  
  84.   /* make sure that n is a power of 2 */
  85.  
  86.   result = fft_binary_logn(n) ;
  87.  
  88.   if (result == -1) 
  89.     {
  90.       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
  91.     } 
  92.   else 
  93.     {
  94.       logn = result ;
  95.     }
  96.  
  97.   /* bit reverse the ordering of input data for decimation in time algorithm */
  98.   
  99.   status = FUNCTION(fft_complex,bitreverse_order) (data, stride, n, logn) ;
  100.  
  101.   /* apply fft recursion */
  102.  
  103.   dual = 1;
  104.  
  105.   for (bit = 0; bit < logn; bit++)
  106.     {
  107.       ATOMIC w_real = 1.0;
  108.       ATOMIC w_imag = 0.0;
  109.  
  110.       const double theta = 2.0 * ((int) sign) * M_PI / (2.0 * (double) dual);
  111.  
  112.       const ATOMIC s = sin (theta);
  113.       const ATOMIC t = sin (theta / 2.0);
  114.       const ATOMIC s2 = 2.0 * t * t;
  115.  
  116.       size_t a, b;
  117.  
  118.       /* a = 0 */
  119.  
  120.       for (b = 0; b < n; b += 2 * dual)
  121.     {
  122.       const size_t i = b ;
  123.       const size_t j = b + dual;
  124.       
  125.       const ATOMIC z1_real = REAL(data,stride,j) ;
  126.       const ATOMIC z1_imag = IMAG(data,stride,j) ;
  127.  
  128.       const ATOMIC wd_real = z1_real ;
  129.       const ATOMIC wd_imag = z1_imag ;
  130.       
  131.       REAL(data,stride,j) = REAL(data,stride,i) - wd_real;
  132.       IMAG(data,stride,j) = IMAG(data,stride,i) - wd_imag;
  133.       REAL(data,stride,i) += wd_real;
  134.       IMAG(data,stride,i) += wd_imag;
  135.     }
  136.       
  137.       /* a = 1 .. (dual-1) */
  138.  
  139.       for (a = 1; a < dual; a++)
  140.     {
  141.  
  142.       /* trignometric recurrence for w-> exp(i theta) w */
  143.  
  144.       {
  145.         const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
  146.         const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
  147.         w_real = tmp_real;
  148.         w_imag = tmp_imag;
  149.       }
  150.  
  151.       for (b = 0; b < n; b += 2 * dual)
  152.         {
  153.           const size_t i = b + a;
  154.           const size_t j = b + a + dual;
  155.  
  156.           const ATOMIC z1_real = REAL(data,stride,j) ;
  157.           const ATOMIC z1_imag = IMAG(data,stride,j) ;
  158.           
  159.           const ATOMIC wd_real = w_real * z1_real - w_imag * z1_imag;
  160.           const ATOMIC wd_imag = w_real * z1_imag + w_imag * z1_real;
  161.  
  162.           REAL(data,stride,j) = REAL(data,stride,i) - wd_real;
  163.           IMAG(data,stride,j) = IMAG(data,stride,i) - wd_imag;
  164.           REAL(data,stride,i) += wd_real;
  165.           IMAG(data,stride,i) += wd_imag;
  166.         }
  167.     }
  168.       dual *= 2;
  169.     }
  170.  
  171.   return 0;
  172.  
  173. }
  174.  
  175.  
  176. int
  177. FUNCTION(gsl_fft_complex,radix2_dif_forward) (TYPE(gsl_complex_packed_array) data, 
  178.                           const size_t stride, 
  179.                           const size_t n)
  180. {
  181.   gsl_fft_direction sign = forward;
  182.   int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
  183.   return status;
  184. }
  185.  
  186. int
  187. FUNCTION(gsl_fft_complex,radix2_dif_backward) (TYPE(gsl_complex_packed_array) data,
  188.                            const size_t stride, 
  189.                            const size_t n)
  190. {
  191.   gsl_fft_direction sign = backward;
  192.   int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
  193.   return status;
  194. }
  195.  
  196. int
  197. FUNCTION(gsl_fft_complex,radix2_dif_inverse) (TYPE(gsl_complex_packed_array) data, 
  198.                           const size_t stride, 
  199.                           const size_t n)
  200. {
  201.   gsl_fft_direction sign = backward;
  202.   int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
  203.  
  204.   if (status)
  205.     {
  206.       return status;
  207.     }
  208.  
  209.   /* normalize inverse fft with 1/n */
  210.  
  211.   {
  212.     const ATOMIC norm = 1.0 / n;
  213.     size_t i;
  214.     for (i = 0; i < n; i++)
  215.       {
  216.     REAL(data,stride,i) *= norm;
  217.     IMAG(data,stride,i) *= norm;
  218.       }
  219.   }
  220.  
  221.   return status;
  222. }
  223.  
  224. int
  225. FUNCTION(gsl_fft_complex,radix2_dif_transform) (TYPE(gsl_complex_packed_array) data, 
  226.                       const size_t stride, 
  227.                       const size_t n,
  228.                       const gsl_fft_direction sign)
  229. {
  230.   int result ;
  231.   size_t dual;
  232.   size_t bit; 
  233.   size_t logn = 0;
  234.   int status;
  235.  
  236.   if (n == 1) /* identity operation */
  237.     {
  238.       return 0 ;
  239.     }
  240.  
  241.   /* make sure that n is a power of 2 */
  242.  
  243.   result = fft_binary_logn(n) ;
  244.  
  245.   if (result == -1) 
  246.     {
  247.       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
  248.     } 
  249.   else 
  250.     {
  251.       logn = result ;
  252.     }
  253.  
  254.   /* apply fft recursion */
  255.  
  256.   dual = n / 2;
  257.  
  258.   for (bit = 0; bit < logn; bit++)
  259.     {
  260.       ATOMIC w_real = 1.0;
  261.       ATOMIC w_imag = 0.0;
  262.  
  263.       const double theta = 2.0 * ((int) sign) * M_PI / ((double) (2 * dual));
  264.  
  265.       const ATOMIC s = sin (theta);
  266.       const ATOMIC t = sin (theta / 2.0);
  267.       const ATOMIC s2 = 2.0 * t * t;
  268.  
  269.       size_t a, b;
  270.  
  271.       for (b = 0; b < dual; b++)
  272.     {
  273.       for (a = 0; a < n; a+= 2 * dual)
  274.         {
  275.           const size_t i = b + a;
  276.           const size_t j = b + a + dual;
  277.           
  278.           const ATOMIC t1_real = REAL(data,stride,i) + REAL(data,stride,j);
  279.           const ATOMIC t1_imag = IMAG(data,stride,i) + IMAG(data,stride,j);
  280.           const ATOMIC t2_real = REAL(data,stride,i) - REAL(data,stride,j);
  281.           const ATOMIC t2_imag = IMAG(data,stride,i) - IMAG(data,stride,j);
  282.  
  283.           REAL(data,stride,i) = t1_real;
  284.           IMAG(data,stride,i) = t1_imag;
  285.           REAL(data,stride,j) = w_real*t2_real - w_imag * t2_imag;
  286.           IMAG(data,stride,j) = w_real*t2_imag + w_imag * t2_real;
  287.         }
  288.  
  289.       /* trignometric recurrence for w-> exp(i theta) w */
  290.  
  291.       {
  292.         const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
  293.         const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
  294.         w_real = tmp_real;
  295.         w_imag = tmp_imag;
  296.       }
  297.     }
  298.       dual /= 2;
  299.     }
  300.  
  301.   /* bit reverse the ordering of output data for decimation in
  302.      frequency algorithm */
  303.   
  304.   status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
  305.  
  306.   return 0;
  307.  
  308. }
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.